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ABSTRACT 


Two different initialization methods were developed and 
tested in global data assimilation experiments covering a 
five-day period. One method was based on the nonlinear 
normal mode initialization, and the other was based on the 
balance equation. Both techniques were developed using the 
calculus of variations methodology. In both methods, the 
initial divergence was computed from the forecast first- 
guess fields, except it was partially modified in the 
nonlinear normal mode method to improve the balance. 

The assimilation system used to test the initialization 
methods was developed for the global forecast model at the 
Fieet Numerical Oceanography Center. This model was adapted 
from the general circulation model developed at the 
University of California at Los Angeles. A comparison of 
the gravity wave noise from the two methods is given for 
versions of the model with and without heating. Other 
comparisons are given for divergence, precipitation rates, 
wave structure and cyclogenesis. [he two methods are 
Similar in their performance in data assimilation. The 
balance equation method is more flexible in weight specifi- 
cation and, consequently the forecasts verify with 


observations closer than the normal mode method. 
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I. INTRODUCTION 


Numerical methods have been useful in the prediction of 
atmospheric flow patterns since the pioneering work of 
Charney et al. (1950) on the ENIAC computer. In that first 
experiment, and in operational models that soon followed, 
pressure height analyses were used as initial conditions for 
the prediction equations. At that time, the numerical 
models excluded gravity waves, and no special processing of 
the initial conditions were necessary. However, when the 
primitive equation (PE) models came into widespread use, the 
independent variables consisted of mass (pressure, tempera- 
ture and geopotential height) and motion (wind, vorticity 
and divergence) variables. Rather than analyze separately 
the wind field, it was computed from the geopotential height 
analyses with a diagnostic field relationship such as the 
balance equation (see Charney, 1955). Since PE models 
permit gravity waves as well as Rossby waves, unbalanced 
initial conditions may become badly distorted during the 
integration. Such imbalances occur when the motion and mass 
Variables are not dynamically matched, as is typically the 
case due to the inaccuracies and the spatial distribution of 
the observations. Therefore, although winds are analyzed 
along with the pressure heights, some method is required to 


remove the gravity wave noise. 


dad 








With the advent of large quantities of asynoptic data 
(data observed at random times, such as from satellites and 
aircraft) came the concept of data assimilation, in which 
the numerical model plays a crucial role of carrying infor- 
mation between observation times (Charney, Halem and 
Jastrow, 1969; Hayden, 1973; Bengtsson and Gustavsson, 1971, 
1972; and Williamson and Kasahara, 1971). The primary 
motivation for data assimilation has been to update the 
numerical model frequently enough that it would represent 
the state of the atmosphere at all times. In this way, the 
asynoptic data could be used in the analysis more effec- 
tively, and thereby improve forecasts. However, several 
difficult problems became apparent from the initial, some- 
what naive, methods of updating a model. The most difficult 
problem has been the inability to assimilate all types of 
data. The mass observations are particularly difficult 
(Kistler and McPherson, 1975; Daley and Puri, 1980; Puri, 
1981) because the geostrophic adjustment mechanism of the 
model tends to disperse unbalanced mass information through 
the gravity waves. Another problem, particularly noticeable 
in the tropics, is that the model may become severely unbal- 
anced. This may set up spurious oscillations on a global 
scale that are very difficult to damp with time filters. 

For example, the pressure field in the tropics may tend to 


Slosh with periods of 12 to 18 hours and amplitudes in 
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excess of 10 mb. It then becomes extremely difficult to get 
the model state to approach asymtotically that of the 
atmosphere. 

Balancing the analyzed data, which removes the gravity 
wave noise, avoids the problems mentioned above. Because 
the wind data may be located in regions without mass data 
and vice versa, the balancing method should be flexible. 
Although much of the gravity wave noise that interferes with 
data assimilation can be controlled with some sort of filter 
applied during the integration, most operational centers 
constrain the initial data so that winds and mass are 
balanced. 

The two basic categories of the balancing procedure are 
the dynamic and the static. The dynamic method involves 
continuously inserting data until the model fields are 
adjusted to the new information. Integration may be 
performed in a forward-backward fashion as discussed by 
Nitta and Hovermale (1969), Temperton (1976), and Haltiner 
and McCollough (1975), or it may be performed in only one 
direction as discussed by Miyakoda et al. (1976) and Hoke 
and Anthes (1976). The main difficulty with the dynamic 
method is inefficiency; it requires the equivalent of a 24- 
hour forecast, and it does not seem to work well for mass 
data (Williamson and Temperton, 1981). Static methods, on 


the other hand, are widely used in operational centers 





(Daley, 1981). In the static method, a diagnostic relation- 
ship is imposed on the analyzed heights and winds. The 
acceptance of the mass data in the model depends on the 
constraints imposed during the balancing (Daley, 1978) or 
whether the analysis of mass variables also produces correc- 
tions to the motion variables (Philips, 1982b). 

Several static initialization methods are available. 
One method utilizes the multivariate optimum interpolation 
to make mass corrections consistent with motion eens then’: 
(Lorenc, 1981; Schlatter, 1975), and then the remaining 
gravity noise is removed with some sort of balancing. 
Unfortunately, this multivariate analysis method links the 
mass and motion through the geostrophic approximation, and 
therefore produces a bias around well-developed systems 
(Williamson, Daley and Schlatter, 1981). Another method 
requires that the mass and motion variables be analyzed 
independently, and then the calculus of variations initiali- 
zation of Sasaki (1958) either adjusts the mass variables to 
the motion variables, or vice versa, depending on the 
expected accuracy in the mass variables relative to the 
winds. The main difficulty with this method is that there 
is no convergence guarantee for the iterative methods 


required to solve these problems (Tribbia, 1981; 1982). 





The constraints imposed on the initial data to remove 
gravity waves are most commonly the nonlinear normal mode 
methods of Machenhauer (1977) and Baer and Tribbia (1977). 
In these methods, the nonlinear component of the balance is 
computed assuming little or no tendency of the gravity mode 
coefficients. No other initialization method is capable of 
suppressing the initial imbalance in a forecast so effec- 
tively. Additionally, the nonlinear normal mode methods 
have the advantages that they provide conditions that are 
compatible with the numerical scheme of the model, generate 
realistic vertical motion in the extratropics, and produce 
balanced flow in the regions with terrain (Daley, 1981). 

In a theoretical study, Leith (1980) used a 
quasigeostrophic model to show that the nonlinear normal 
mode methods are nearly the same as constraining the Hite el 
conditions to the balance and omega equations. Therefore, 
it seems reasonable to expect that the balance equation 
might still be competitive with the normal mode methods, 
particularly since the balance equation constraint is 
relatively easy to impose. 

There remains the problem of generating an appropriate 
divergence to go with the nondivergent winds produced by the 
balance equation. Tarbell et al. (1981) used a modified 
omega equation which improved the precipitation forecasts of 


a mesoscale model during the initial hours. Considering 


IPS, 





that the omega equation, and even the nonlinear normal mode 
method do not generate divergence patterns in the tropics 
that are compatible with the latent heat release (Bengtsson, 
1981), the divergence from the forecast first-guess may be 
the best estimate for divergence. 

In this study, we examine the adequacy of the forecast 
first-guess divergence. Because this PSGene is gene- 
rated by the model, it is compatible. Thus, the interruption 
to the cumulus convection in the tropics is minimized during 
assimilation. Another goal of this study is to determine 
whether the classical balance equation method of Charney 
(1955) used in the variational method of Sasaki (1958) is 
competitive with the more elegant, yet cumbersome, nonlinear 
normal mode method. Since one of the main benefits of the 
nonlinear normal mode method is the divergence it generates, 
comparison of the balance equation and normal mode method 
is, in many respects, a test of the accuracy of the forecast 
first-guess divergence. 

In the following, a brief description of the data 
assimilation system used in this study is given. Besides 
the initialization methods, this system includes an 
objective analysis method for wind and pressure height 
observations and the global finite difference model 
described by Arakawa and Lamb (1977). The initialization 


methods, which are the primary topics of this report, are 
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presented in detail. Both methods are presented in a 
calculus of variations framework. The normal mode method 
uses Machenhauer's (1977) initialization as a constraint, 
whereas the other method uses the balance equation as a 
constraint. Finally, several data assimilation experiments 
are described that illustrate some of the characteristics of 
the different initialization methods, particularly with 
regard to precipitation rates and divergence during the 


early forecast hours. 
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IT. THE DATA ASSIMILATION SYSTEM 


In the data assimilation method used in this study, the 
prediction model is periodically updated at 12-h intervals. 
Each update requires several steps. First, the 12-h fore- 
cast is interpolated to the analysis coordinates, which are 
Prescmme SUrtaces on.a 2.59 by 2.59 grid. This forecast 
becomes the first-guess field. The objective analyses of 
wind and pressure height are done with a three-dimensiona|l 
successive corrections method (See Appendix A) on the 
Standard pressure surfaces (50, 100, 150, 200, 250, 300, 
400, 500, 700, 850 and 925 mb). However, the surface 
pressure analysis is produced from the calculus of varia- 
tions method of Hol] and Mendenhall (1972). At this time, 
the initialization is conducted. The balance equation 
initialization is done on pressure surfaces prior to the 
interpolation to model coordinates, whereas the normal mode 
initialization is done in model coordinates. These initial- 
ization methods are described in the next two chapters. 
Finally, a 12-h forecast is made in preparation for the next 
update. 

The other variables in the model, such as boundary 
layer depth and moisture, are not updated with new data. 
Rather, they are carried forward from the forecast first 


guess. 
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The numerical model used to assimilate the data was 
developed at the University of California at Los Angeles 
(UCLA) to study the general circulation of the atmosphere 
(Arakawa and Mintz, 1975; Arakawa and Lamb, 1977). In this 
model, the resolution of the mass variables (surface pressure 
and temperature) is 4° latitude by 59 longitude and six 
levels from 50 mb to the surface. The zonal wind is defined 
one-half grid interval to the east and the meridional wind 
is defined one-half grid interval to the south (Scheme C in 
Arakawa and Mintz, 1975). The heating parameterization 
includes the Arakawa and Schubert (1974) parameterization 
interacting with a bulk parameter boundary layer (Randall, 
Be7ij; Lord, 1978). 

The time differencing is a combination of five leap- 
frog steps for each Matsuno backward step, while the heating 
is computed during a single forward step preceding each 
Matsuno step. 

In the next section, the model's normal modes are 
computed, and more detail is given on the dynamical forecast 


scheme. 
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IIT. NONLINEAR NORMAL MODE INITIALIZATION 


Before nonlinear normal mode initialization can be 
applied, the normal modes of the linearized model equations 
are required. These are obtained below by separation of 
variables of the linearized equations. Unfortunately, the 
normal mode methods match the motion variables more closely 
than the mass variables (Daley, 1981). To establish control 
over this mass rejection mechanism, the normal mode method is 
converted to a variational one similar to that of Daley 


(1978) and Tribbia (1982). 


A. VERTICAL MODES 
The linearized governing equations may be written with 


the vertical component in vector form as 


f) - m 
s¢ V fee ke Wt VCRT In Bee t ¢) = Q WV? (Salk) 
fe) = 
mee ee Oo ae) 
2 inp.+ mg. vV) = Q., and (ene 
at Sk : pe 
> ee : U 


Here W is the vector form of wind,T is temperature and IT is 
the rest state temperature. Pe. is surface pressure, 9» iS 


geopotential, v-Wis divergence and 9. is terrain 
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geopotential. +t, 2 and G are linearized operators defined 


~~ 
~~ 


in Appendix 8B. Vy QT and Qp are the nonlinear components 


of their respective equations. 

Following Temperton and Williamson (1981), the 
temperature and surface pressure may be described by a 
Single variable by using 

gh = + RT In D ; (325) 
Operating on (3.2) with G, multiplying (3.3) by RT, and then 
adding the resulting two equations gives a single equation 


for mass, 
g bh +G c(v: W) + RTC (v2 W) = GQ, + Q,RT (3.6) 


which may be rewritten as 


where 
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Time Cama im set (3.7) is vertically coupled, but by 


separation of variables, it can be transformed into a set 


that is not coupled. This is done through the identity 
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where matrix — and diagonal matrix D contain the eigen- 


vectors and eigenvalues of C, respectively. Transforming h 


md W, or 
W= 87) and Mure 
= eo h , (3.12) 


produces the following uncoupled equation set: 


i" 


W + fikevKW + 9 (h) = Qy and (3.13) 
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ct 
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where y and Qh are transforms of Qy atid Qe, “respectifely. 

The independent variables in (3.13) _ (3.14) are the 
coefficients of the vertical modes, i.e., the eigenvectors 
c. Naturally, there are as many modes as there are levels 
i the model. 

The eigenvectors, &, are shown in Figures 1 and 2 for T 
equal to a constant (3009X) and for T equal to (209, 214, 
233, 254, 270, 281)9K. The corresponding eigenvalues 
(equivalent depths) are given in Table 1. Notice that while 
the profiles do not significantly change shape for the 
different temperature profiles, their eigenvalues are 
considerably different. Also, the vertical modes of tne 
Arakawa and Lamb (1977) model shown here are noticeably 
different from those given by Temperton and Williamson 


(1981) for the same values of T. Specifically, the peaks in 
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Table 1. Equivalent depths (eigenvalues) of the vertical 
modes of the Arakawa and Lamb (1977) model. 


Level] T Constant T Average of 
300°K 1 March 1965 
FE 9,291 o 7,874 m 
2 1,689 784 
3 387 E71 
4 126 55 
5 45 18 
6 10 c 


the profiles near the model top are not present in the 
gravest modes given in Figures 1 and 2. Since these modes 
are insensitive to the values of rest state temperature, 
they need not be updated for each new data set. 

The equivalent depths are sensitive to the location of 
the model top. For example, changing the model top from 50 
mb to 0 mb increases the equivalent depths of the externa! 
mode from 7/874 m to 9660 m. The consequence of keeping the 
model top at 50 mb is that all the vertical eigenvalues or 
equivalent depths are smaller than if the top was placed at, 
say, 10 or 20 mb. This becomes a factor when the nonlinear 
balance is performed, as will be discussed later. 

After h is balanced, the inverse of (3.5) is required 
ta solve for re and T. This is done following the approach 
of Temperton and Williamson (1981) and Andersen (1977), 


where y-y may be eliminated between (3.3) and (3.7) to give 
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G2 yg = ONET. (9.45) 


The nonlinear components are grouped in the term NLT. From 
linear theory, it is possible to relate h adjustments to 
Ps adjustments, or 

A in ) a al (Be q Any, (3.16) 
where the adjustment balances the analyzed fields. 

In a similar manner, elimination of V-W between (3.2) 
and (3.7) produces 

AT = 2 con") g Ah Cay) 
Using definition (3.5), it is easily shown that (3.16) and 
(3.17) are consistent methods for determining Alnp, and aT 
from Ah. 

A two-grid interval wave was found to exist in the 
solution of (3.17). To eliminate this problem, a varia- 
tional method was used which minimizes AT, thereby removing 
the two-grid length waves. In this method, Ag is computed 
from (3.5) and then AT is determined from the minimization 
of 


(AT ie + B(A i + 2r-(Ad, - ‘ 
a ee Teele 2 


io ae 


GAT, ) Ag, - (32 13) 


Le 1 


The details of this solution are found in Barker (198lc). 
This method reduced the size of the root mean square (RMS) 


corrections to about one half of those obtained using 


hay). 





In this section, the vertical coupling of the 
linearized equations (3.1)-(3.4) was eliminated through 
separation of variables. Following a similar procedure, the 
horizontal coupling can also be removed, as shown in the 


next section. 


Br. HORIZONTAL MODES 
The equations (3.13) and (3.14) in linearized finite 


difference form for each vertical mode are 
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The mode index, m, is assumed for D and each variable. 
The finite difference operators are 

C8 T)y = Tke1s2 + Tko1/2> 

me ieee * 1 ee1/2 
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The other variables are defined as follows: a is the earth's 
radius, Ad is the longitudinal grid interval, Ae is latitu- 
dinal grid interval, u and v are the east and north compo- 
nents of wind, respectively, i and j are the longitudinal 
and latitudinal indexes, respectively, and o is coso. 

The linearized equations (3.19)-(3.21) are derived from 
the model equations in the flux form given by Arakawa and 
Lamb (1977). However, linearization of these equations 
removes their flux form, and the only terms absolutely 
unique to the model are the Coriolis terms, which differ 
from other models in the way that f is averaged. As it 
turns out, special definitions for the Coriolis terms are 
desirable in order to keep the matrix operator of (3.19)- 
(3.21) symmetric, which simplifies the corresponding 
eigenvector matrix so that it can be inverted with a 
transpose operation. This decreases the computer storage 
and time requirements needed to transpose the variables into 
normal mode coefficients to one half that required for non- 
symmetric operators. Fortuitously, the small errors 
introduced by these modifications do not affect the normal 
mode balancing (Temperton and Williamson, 1979). To achieve 


symmetry, the Coriolis term in (3.19) is replaced by 
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— — te 7 
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and the Coriolis term in (3.20) is replaced by 


Pe ag i ae 
- (u),; . + fey Cu) = 
2 
where 
et OUEST eae ee, 
J cos (“3 
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= 1/3 Fieis2 + 2/3 f5_1 72 
J cos (FH) 


These definitions correspond to a potential enstrophy 
conserving finite difference scheme as derived by Temperton 
and Williamson (1979). 

From this point, the procedure closely follows that of 
Dickenson and Williamson (1972) except that the finite 
differences are written for scheme C (see Arakawa and Lamb, 
1977; Temperton and Williamson, 1981). The dynamical state 
vector is defined as 


u(y. 7a 52M) 
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By assuming a wave solution in the form 
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and consequently, 


™ Le 
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(3.19)-(3.21) become 
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where 
m'(k) = cos(Ss)-i sin(S+), 
Ke OS sin( SSA) /(4A) , r(k) = cos (<4) , and S ¢ 1s the 


filtering applied near the poles to keep the gravity-wave 
terms computationally stable (Arakawa and Lamb, 1977). 


The symmetric operator 
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or in Watrix form, 
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For the general case when j is not near the north pole or 


the equator 
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To produce continuous Coriolis forcing near the poles, 
a computational u is defined. At the North Pole, 


Luly = 0, where 


the [ ] represents a zonal average. 


Continuity requires that 
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The continuity equation at the pole is 
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A grid with a four degree meridional interval puts the 
equator at a v-point. Using Y1/2 as the equator point, the 
equations can be written for the symmetric and antisymmetric 
components. For the symmetric component, U, = Ups hy = No 


and V1 /2 = 0. The h- and u-equations are 
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Changing the model resolution will generally require that 
these equations also be changed, depending on how the 
variables are staggered relative to the equator. 

To complete the computation of the horizontal modes, 


the matrix equation (3.32) is rescaled using 
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This is the wave equation of the expansion Coertriciemes, C, 
that we have been seeking. The elements of C are functions 
of the vertical, zonal and meridional mode numbers m, k and 
1, respectively. These coefficients are the amplitudes of 


the various modes required to represent a particular atmos- 


pieric state, i.e., 
. (3.46) 


The nonlinear term is now r, and the mode frequencies are A. 
For each m and k, there are 3N equations for C; 2N are 
gravity waves and N are rotational waves. N, in this case, 
is the number of degrees of freedom in the meridional 
direction. 

The structure of the modes from (3.43) is given in 
Figures 3 and 4. They are nearly identical to the modes 
published by Dickenson and Williamson (1972). The scaling 
is different, but this is of no consequence, as the modes 
are put in orthonormal form before they are used. The 
frequencies of the various modes corresponding to those 
given by Dickenson and Williamson (1972) and Temperton and 
Williamson (1981) are given in Tables 2 and 3. Resolution 
differences account for most of the variability of the 
frequencies, which can be seen from the computations of 
Dickenson and Williamson (1972) for two different resolu- 


tions. Their 5° unstaggered grid results are similar to the 
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Fig. 4. Similar to Fig. 3 except for selected 
gravitational modes. 
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maole 2. 


10 
Nh 
Na 
8 
14 
1s 
16 


Frequencies of Rossby modes for D = 10 km, k=l 
from computations of Temperton and Williamson 

(1981) (T&W), Dickenson and Williamson (1972) 

(D&W), and for the model used in this study, B. 
The horizontal grid intervals are specified in 
degrees. 


T&W, 10° D&W, 5° D&W, 2.59 Bee 

6.11 E-05 6.12 E-05 6.14 E-05 6.14 E£-05 
1.44 £-05 1.43 £-05 1.45 £-05 1.45 £-05 
8.64 E-06 8.60 E-06 8.73 E-06 8.74 £-06 
5.72 E-06 5.73 £-06 5.87 £-06 5.88 £-06 
3.98 E£-06 4.02 E-06 4.17 E-06 A eG 
2.87 E-06 2.93 E-06 3.08 E£-06 3.09 £-06 
2.14 E-06 2.20 £-06 2.36 E£-06 2.36 E-06 
1.63 £-06 1.69 E-06 1.86 E-06 1.86 £-06 
1.27 E-06 1.32 E-06 1.49 £-06 1.49 £-06 
1.01 £-06 1.04 £-06 1.22 £-06 1.22 £-06 
8.10 E-07 3.18 E£-06 1.02 £-06 1.02 £-06 
6.62 E-07 6.43 £-06 8.58 E-07 Soy ee =O 
5.52 E-07 4.99 E-07 7.30 E-07 7.30 E-07 
4.70 E-07 oa E207 eal Say 6.28 E-07 
4,11 E-07 2.71 £-07 5.43 £-07 5.45 £-07 
3.75 £-07 1.75 E-07 4.73 E-07 4.76 E-07 
3.13 £-07 8.60 E£-08 AA E07 4.18 E-07 
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~5.44 
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cs esy) 
S75 
=2\,78 
2222 
-3.63 
-4.01 
-4.36 
-4.69 
-4.98 
aee25 
~5.44 
-5.61 
Se? 
~5.94 
-5.94 


Similar to Table 2 except for frequencies of 


eastward gravity modes for D 
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E-05 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
E-04 
F-04 
E-04 
E-04 
E-04 
F-04 


D&W, 
3.6 
0 
~ 85 
ors 
7.6 
so) 
61 
. 00 
~ 34 
A657 
295 
Rall 
~42 
ao 
.68 
FD 
mon 


50 
E-05 
E-04 
F-04 
F-04 
E-04 
E-04 
E-04 
F-04 
E-04 
E-04 
E-04 
E-04 
F-04 
F-04 
F-04 
E-04 
E-04 
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= 10 km, 

a 

E-05 -5 
E-04 -1 
E-04 -1 
E-04 -2 
E-04 -2 
E-04 -3 
E-04 -3 
E-04 -4 
E-04 ~4 
E-04 -5 
E-04 -5 
E-04 -6 
E-04 ~6 
E-04 -6 
E-04 -7 
E-04 -7 
E-04 -8 


kK=1. 

By 4°x5° 
e338SE=00 
73S07E=08 
~87 E-04 
.36 E-04 
~84 E£-04 
~31 £-04 
7/8 E-04 
.25 E-04 
7/1 E-04 
~18 E-04 
~.64 E-04 
.09 £-04 
~54 E-04 
.99 £-04 
-43 E-04 
-86 E-04 
-29 E-04 





Temperton and Williamson (1981) 10° staggered grid values, 
whereas their 2.59 nonstaggered grid results are similar to 
those computed from the above equations for a 49X59 
Staggered grid. As Temperton and Williamson (1981) point 
out, the staggered grid produces modes comparable to those 


of a nonstaggered grid with half the resolution. 


Cc. NONLINEAR BALANCING 

Machenhauer (1977) discovered that the nonlinear terms 
have a much slower time variation than their respective 
gravity modes. Under this first order approximation, the 


equation for the gravity modes from (3.45), 
se C(k,2,m) = -ivC(k,2,m) + r(k,2,m) (3.47) 


has the solution (3.48) 


) 9 3 =] t 
C(k,2,m,t) = a. + Caeaaamms )) = rik tom) ad 


Therefore, removal of the fast modes requires that the 


second term on the right-hand side of (3.48) be zero, i.e., 
Cock ,2 amo) = r(k,2,m)/(iv). (3.49) 


The subscript B signifies the balance condition. 
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Rather than define the nonlinear terms, it is easier to 
use the model to make this computation. This is done by 
running the model for one time step and determining the 
tendency of the gravity mode coefficients. The nonlinear 


term is the total minus the linear tendency, or 


mik, Ua) = Clk t,matj-C(k,2,m,0) t iv C(k,2,m,t). 
(37-50) 


Ve is the computational frequency for the forward timestep, 


which may be determined from the analytic frequency by 


E 1 2 
Ve = arctan Nae ae Tar bg Meet a (Corsa) 
using standard methods. Using v,. for v in (3.49) and then 
substituting (3.50) gives the corrections to the fast modes 


required to balance, or 
ACR = -i(C(k,2,m,at)-C(k,2,m,0))/ (atv (k,2,m) ) (395'<) 


Leith (1981) used the quasi-geostrophic relations to show 
that this balance is equivalent to using a slightly modified 
version of the balance equation along with the omega equa- 
tion. Phillips (1981) points out, however, that more than 
one iteration with (3.52) introduces new terms not 


consistent with quasi-geostrophic balancing, and therefore, 
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should be avoided. The Baer-Tribbia (1977) method, however, 
does allow more iterations, and will be tested in future 
versions of the normal mode balance. 

Additionally, Machenhauer's method (3.52) does not work 
except for those gravity modes with periods equal to or 
less than any of the rotational modes (Ballish, 1981; 
Phillips, 1981). For the version of model used here, this 
restriction limited the modes that could be balanced to all 
of the external and first internal, and afew of the second 


internal, modes. 


vy. VARIATIONAL BALANCING 

One of the difficulties with the nonlinear normal mode 
method is that the rotational modes are primarily determined 
from the vorticity of the analyzed wind fields (Daley, 1981; 
Daley and Puri, 1981). This becomes a problem over vast 
regions of the globe in which the principal observations are 
remotely sensed temperature profiles from satellites. 
Therefore, to describe these regions adequately, the 
initialization must be able to assimilate mass observations. 
A possible method for incorporating mass observations is to 
convert the normal mode balance into a variational framework 
(Daley, 1978). Tribbia's (1982) spectral shallow water 
approach was adapted for this purpose, except nere the 
method is developed for a multilevel finite difference 


model. 
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The normal modes introduced by (3.43) may be written 


exp(-ika,)/I Ve (Zee) 





but instead of (3.46), the coefficients of the rotational 


modes are determined from a specially designed inner product 


C(k,2) =Cy-¥"(k, 2)> (3.54) 


where 


Crete) 


( 3a055 ) 


I J z 
re < (6, Ag aq sg) Cues rpaq/2)U(2.9, 8.)a (8 ) 


st 


War s2rtq ie Vos oe ) V(2,0, j) doles 472) 


+. 


We (0505 )Lo( 9502, )o(2, 8; )]o(e ,)bexp(-ika,)/1'/* 


Wy and We are the weights for winds and mass, respectively. 
o is the same as in (3.25)-(3.27). Because no attempt is 
made to put a vertical variation in the weights, the 
vertical mode index has been dropped from Teas Cone Tere 


The inverse of (3.46) is 


(3.56) 
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If yy is the balanced initial condition and Yo the 
analysis, then an optimal initial condition is that 


obtained from a minimization of 


L = (yy - yg)ely - 9) (3.57) 


~ 


with respect to the modal coefficients. 
The balanced data have components on both the rota- 


tional and gravitational manifolds, i.e., 


ie (ne Yc: (3.58) 
where 
Yo = 2 ae) f and 
I-] G 
” K <k 
= © Ff Gry 

: 3.60 
~2 k=O en] 2% -* | | 


xk and Gk represent the rotational and gravitational mode 
coefficients, respectively. Notice that G plus R equals N, 
which is the total number of modes. The nonlinear balance 
relationship requires that the gravitational component be a 
function of the rotational component only (Phillips, 1981; 
Baer and Tribbia, 1977; Daley, 1981), i.e., GX = G(X), 


therefore minimizing (3.57) requires that 


ce = (3.61) 


I 
7X8 
04 


5 | 





for the rotational mode coefficients X, 


Ms Tribbia (1982) notes, 


TMi.S 
this 


poss 


The 


thas 
1S ert ie harly tnue for 
equation more tractable, 


ible where the first pass 


ae 


(1) 
I 


t-<> 


Superscript is 


where 


ites ation» cycle. 


or 


(3962) 
equation is not easily solved. 
non-zero values of G. To make 


an iterative solution is 


is solved using 


(3.93) 


This sigplifies (3862) 


qe Ext et incity-e oe 
reduces to a linear equation set 
A Dts 7 (3.65) 
AX = an: “ ee and 
~ kK=0 2=1 a : 
; <7 


for 


al] le Gmecthne rotataania lem 


anifold. 








If the weights are unity everywhere, then the orthonormality 
of the modes makes A an identity matrix, and YI is simply 
equal to the rotational component of We 

subsequent values for X may be computed using (3.62). 

G is evaluated from Machenhauer's equation (3.52) using the 
value of X from the previous cycle. The variation of G with 
respect to Ke is computed numerically. The benefits of 
following this complex procedure are almost certainly not 
worth the effort. 

To avoid the computational workload of (3.62), 
the variational balance is performed on the analyzed correc- 
tions. Assuming the corrections are much smaller than the 
total analysis, the approximate equation (3.64) is more 
accurate than it would be if the balance were performed on 
the total analysis. 

The solution is still not very easy. For total 
variability of the weights, the equation set (3.65) 
represents (I/2°R) or 1656 equations. The corresponding 
Biize of A is over 2.5 million elements. M[herefore, the 
computation of the inverse of A is extremely cumbersome, 
especially when the computer memory available is only about 
100,000. Restricting the variations of the weight to 
latitude only decreases the size of A to 23X23. But now we 
must solve 216 different arrays, one for each vertical and 


horizontal mode configuration. These are solved only once 


and stored. 
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Some attempts were made to add longitudinal variation 
to the weights by iterating, but this did not work. 

Firstly, it was not possible to keep Y] on the rotational 
manifold. Secondly, the variations in the magnitude of the 
weight required to influence significantly the result were 
too large for convergence. 

Ultimately, it is hoped that a subset of the rotational 
manifold may be found that still defines the important 
meteorological information that needs to be retained by the 
initialization and yet has reasonable array sizes. 

A summary of how the final method works is given bya 
Leith (1981) manifold schematic in Fig. 5. Prior to the 
update, the model is assumed to be on the slow manifold (M), 
Say at point 0. Adding the rotational component puts the 
model on the data manifold. This step is represented by the 
line between O and 1, and is equivalent to linearly 
balancing the corrections and inserting them into the model. 
Note that only the rotational modes are affected by this 
step, where (3.65) is solved for X and then the rotational 
component of the corrections is determined from (3.59). The 
imbalances introduced by updating are removed by solving 
(3.52) and then (3.60). This is the nonlinear step and is 


represented by the iine between 1 and 2. 
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Oh, 
a, . 


Gir. 5. Manifold schematic of the update and balance 
Procedure used in the data assimilation. 


Point 0 
represents the forecast before the update. 


The 
path between 0 and 1 represents the linear balance 


step. The path between 1 and 2 represents the 
nonlinear balance step. 
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EB. TESTS OF THE NORMAL MODE METHOD 

The magnitude of the rotational and gravitational com- 
ponents of the analyzed corrections were computed from real 
data taken from the 00 GMT 16 Nov 1979 analysis (see Figs. 6 
and 7). The surprising thing about these curves is that the 
gravitational component is as large or larger than the 
rotational component. This is especially true of the 
Surface pressure. In fact, the surface pressure gravita- 
tional component is so large that the two components must be 
of opposite signs over much of the globe. Because the 
nonlinear component of the initialized fields is determined 
from the rotational component only, the gravitational 
component of the analysis is not used and therefore repre- 
sents lost information. It is also a measure of how well 
the analysis dynamically matches the wind to the mass. 

Completeness theorems for the normal modes require that 
the sum of the rotational and gravitational components of 
the data be equal to the original data. As a test of the 
approximations and final code, the degree to which the 
completeness theorem applied was determined. This was done 
by computing the rotational component of the corrections, 
and then computing the gravitational component from what 
remained of the corrections after the rotational component 
was removed. The corrections were then reconstructed by 


adding the rotational and gravitational components. These 


56 





"€ [LAAaL Sayeubtsap ¢€ Ydtuosqns ‘(aul| 





paysep) yuauodwod [Leuotzeqyou Sszt pue (aul{ pttos) susAjeue a1yepdn (oh S15 
apnz13e) opnd}Ie] 
0°O6 8°e9 0°OL Q°G 0°eL- @°09- @°As- 0 °G6 8°09 8°Of 6°e @°0£- @°09- 8 °O6- 
e0°O 60°8 
ie a ‘ oe 
a“ Ps iis “" ‘ a Ta 
. a7. ao ’ ‘Vat ‘ Laer ~— 
a ard Unt ye > 
— 
Come 
Q0°S 08°S 
epnjz!ie7 opnziyey 
0°86 6°89 @°O£ 8° @°OL- @°69- @°06- 8°e6 @°e9 @°C£ 9° Q°Rz- @°as- 0°0R6- 
a ae ea Tet o°9 
A 
= 
Ta) 
> 
< 
wo 





ae 





TO0) (OUE2 ](s)) iG: BPsstoasqe ayy ueau 

oult Ppoysep yuep ayy Aq uMOUS SL plosLUeW [PuOLze JOU UO [RUOLZeYLAPUH 
OY} AIYZLI OFUO PaALOSaU You BuNnssaud JICPJUNS pue sunzyeuaduay 

94} JO LENPLSAu ayy *Lapow ayy jo doy AU} WOUJ LAADL Pulyy suog 





(9UL{ paysep) syuauodwod [PuOLJePLAeCUH “Lay. pue sSuotidd4u0D pazA,euy +) 2s 
117 
Fanstivy7 ani | 
; ac. "AO-~ "A6- : ° ° a 0°@L£- @°O9- 8°86- 
@°86 e-es @°O£ e"0 Q°O£ 0'e9 @ —_ 0°06 0°e9 e°0f or 
a > 
no 4 
A 
= aa: 
Lo°0—— 0's n0°S & 
Eh’ --- = Oo 
BO o 
JONLILVI JGNLILVI 
0°06 8°e9 @'O£ ee @°@E- = @°09- 9°@6- 0°06 e'e9 6°02 e°@ oC. ee o's 
e°0 
> 
Cc 
> 
% °] 
a 
Mand * 
e°O! = e°et : 
» 
= 
a 0 
@ 
(@) 





98 








reconstructed data werecompared to the original data by 
computing RMS differences (see Fig. 7). The differences 
between the two wind variables were essentially zero, but 
nonlinear relationships between surface pressure and temper- 
ature caused significant residuals to remain in the mass 
variables. The long dark dashed lines near the abscissa in 
Fig. 7 (c) and (d) show the residual for level three of the 
model and the numbers above the curves give the global RMS 
averages. It can be seen that the temperature residual is 
largest, and represents about 14% of the total correction. 
The pressure residual is only about 2%, however. 

As a test of the variational technique, the linear 
balance was rerun with different weights on the mass varia- 
ble each time. Figure 8 shows the global RMS residuals 
between the balanced and analyzed variables for the differ- 
ent weight values. Clearly, the most sensitive region is 
for weight values between Q andl. For aloss in the wind 
residual of about 0.4 Mm sect, the improvement in the mass 
variables is 0.39C and 1 mb. Increasing the weight from 1 
to 10 does not dramatically improve the fit of the balanced 
mass variables to the analysis. Eventually, for very large 
weight on geopotential, the wind residual becomes as poor as 
the forecast first guess. I[t is interesting to note the 
limit in the residual of the wind variables. Even when the 


weight on the mass variables is zero, the wind residual is 


5, 





F deg). 
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RMS( AP} 
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The RMS differences between the analyzed and 
initialized variables for different weighting on 
the mass analysis, Wg. The averages are for the 
entire globe and all levels. The dashed straight 
line shows the RMS differences between the 
analysis and the forecast first-guess. 
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2.8 m sec~t. The reason for this is that the analyzed wind 
contains unrealistic divergence that is orthogonal to the 
rotational manifold (Daley, 1980), and regardless of the 
weighting, this component has no influence on the rotational 
component. 

As a further demonstration of the variational method, 
the same analysis was balanced two different ways. In one, 
the mass variables were weighed very heavily (Wd = 1000) 
poleward of +309 latitude. In the other, they were given no 
weight. The 500 mb vorticity and geopotential fields are 
Shown in Figs. 9 and 10 with the analyzed contours for these 
two cases. From these figures, it is clear that this range 
of weighting is sufficient to map exactly either vorticity 
or geopotential onto the rotational manifold. This is 
probably not a good method for computing mass from motion 
and vice versa, because the relationships are linear. In 
fact, winds computed from the method compare less favorably 
with observations than the 12-hour forecast first guess. 

To correct this problem, Phillips (1982a) proposes a method 
whereby the nonlinear component is removed from the analysis 
prior to balancing. 

The variational method described in this section is too 
cumbersome to impose full variability into the weights. As 
a result, restrictions that only allow variations with 


latitude were imposed. Such limitations are only necessary 
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oh Updated 500 mb vorticity when (a) weight on mass 
is zero, (b) weight on mass 1s 1000 poleward of 
30°. The unprocessed analysis is dashed and the 
contour interval is 25-1076 sec7?t 
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ee ae ee ! el ee — 
Updated 500 mb geopotential when (a) weight on 
mass is zero, (b) weight on mass is 1000 poleward 


of 309. The unprocessed analysis 1s dashed 
and the contour interval is 60 m. 
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in the normal mode method. As shown in the next section, 
the variational method that uses the balance equation as a 


constraint is capable of using weights with large spatial 


variation. 
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Iie THe eBABANCE EQUATION INITIALIZATION 


.. THE VARIATIONAL FORMULATIONS 

For simple models, the nonlinear normal mode balance 
is nearly equivalent to applying the balance and omega 
equations as constraints to the initial conditions (Leith, 
1980). In the variational method developed by Sasaki (1958, 
1970), Stephens (1972) and Haltiner et al. (1976), and 
described in this section, the balance equation is utilized. 
But, instead of the omega equation, the vertical motion is 
derived from the forecast first-guess divergence. 

To impose the balance equation as a constraint, the 
functional 

I($,¥) = Siie=s)? + 3( W- W,)* + 2,g m(u,o) dA (4.1) 
is minimized. Here, » is geopotential, W is the vector 
wind, vw is stream function, A is the horizontal area over 
which the integral is applied, and Ag is the Lagrange 
multiplier. The constraint is 

Myo) = v-fuy + 2d(u,v) - 7%» = 0, (4.2) 
where u and v are the east and north components of wind, 
respectively, J is the Jacobian operator and all operators 
are applied in spherical coordinates. The minimization is 
achieved when the first variation of I{o,u) with respect to 
d5Ap and vw is set to zero. Neglecting the variation of 


J(u,v) and noting that 


Tala A =) 20 (43) 
A 
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for integration over the globe gives 


a 2 
tom Yoke (4.4) 


and 


Z 


BV“py = veB-V(y-y) +gy° 


v + yf.vapg t+ fv“n3. (4.5) 


The Euler-Lagrange equations are (4.2), (4.4) and (4.5) and 
the unknowns are »y, » and yg. The solution of the Euler- 
Lagrange equations is accomplished by eliminating y and » in 


(4.2) using (4.4) and (4.5). This gives: 


2 
Ti (arg)” = E o@(argy™ = mort yt, (4.6) 
(ao)” = Pia” and (4.7) 
peed n n 
f + Vf- 
Mirage. 1 Akg) 7 PFC Ada) (4.8) 


B 


The superscript, n, is the iteration count and the delta 
Specifies the difference between the current and previous 
ioeiation, 1.e., 

Go) So See (4.9) 
where 

Nea = b when n=1. (4.10) 

Note that all boundary conditions are automatically 

eliminated when the solution is desired for a full sphere. 
The second term on the right of (4.8) was dropped because it 
caused unrealistic adjustments to the winds near the 


equator. Consequently, only elliptic equation solutions for 
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7(aa) and Aw are necessary to solve (456)-(4.8). Each 
iteration through (4.6)-(4.8) reduces the residual of the 
balance equation constraint by an order of magnitude, 
therefore two cycles produce sufficiently balanced results. 
A direct solution technique for elliptic equations was used 
to solve (4.6) and (4.8) (see Swarztrauber and Sweet, 1975: 
Rosmond and Faulkner, 1976). 

The forecast first-guess divergence is retained at the 
update time in the following manner: First, the variables 
are balanced as described above. Secondly, the rotational 
component of the forecast first-guess wind field is computed 
and subtracted from the balanced wind. Finally, the 
Variables are interpolated to model coordinates where the 
wind, which is now a nondivergent correction, is added to 
the forecast first-guess wind field. The mass variables, 
however, are balanced and interpolated to model coordinates 
in the conventional manner. In this way, the erroneous 
divergence near steep terrain regions remains small, and the 
first-guess divergence is untouched. 

Another method of retaining the forecast first-guess 
anvergence 1s to balance the corrections rather than the 
updated values. This makes the nonlinear term a little more 
cumbersome, but this procedure has the advantage of not 
affecting areas without new data. In this case, the 
nonlinear term is defined as 


eee ee ey Uy) (4.11) 





where the prime is used to designate a correction value 
rather than an updated value. The integral to be minimized 
TS 

I(g'sy') = Soy? + a( W')© + 2x8 m(w'so') GA, 
where the constraint is 

m(w',o') = vefyw' + 20'(u,v) - 7~6' = 0. 
This basic approach has also been proposed by Phillips 


(1977). 


Br, VERTICAL FILTERING WITH EMPIRICAL ORTHOGONAL FUNCTIONS 

The major problem with the foregoing formulations is 
that the temperature adjustments are not small everywhere. 
For example, if balancing caused the 925 mb height to 
Change 10 m and the 1000 mb height did not change, the 
resulting temperature correction between these two levels 
would be 49°C. Such inconsistencies are difficult to avoid 
for grids covering very large regions. 

To minimize the effects of these inconsistent modifica- 
tions in the vertical, the variables are vertically coupled 
by projecting them onto basis functions prior to the 
balance. Because filtering requires that some of the 
unimportant components be dropped, the empirical orthogonal 
functions are the most suitable functions for this purpose. 


Tnese functions are derived so that they form abasis set 
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that produces the least error when a partial series is used 
to describe a particular three-dimensional field, such as 
geopotential (see Appendix C). 

An example of the efficiency of the empirical 
orthogonal functions is shown in Fig. 11 for a ten-level 
analysis of geopotential. Here, the first mode explains 95% 
of the total variance, whereas the first two modes explain 
98%. Only four modes are necessary to explain 99.9% of the 
total variance. Projecting the wind analysis onto the first 
four modes retained the details of the jet streams produced 
by the analysis. Holmstrom (1963) first noticed this very 
rapid convergence of the empirical orthogonal functions in 
representing geopotential profiles. The smoothness of the 
first four empirical orthogonal functions insures that 
inconsistent vertical variations of wind or geopotential 
between adjacent levels will be truncated when used to form 
a partial series of the analyzed fields prior to balancing. 
Another advantage of this method is that considerable 
computer time is saved. Rather than treating 10 levels, 
only the four vertical modes need to be balanced. 

To show how empirical orthogonal functions are 
incorporated into the variational balance, the balance 
equation is written in vertical vector form 


mv.) = V-(fvy) + A - Veo, ala) 


69 





SECOND 





2.9% 


FOURTH 





0.017% 0.005% 


ree + 


Eudes site First six empirical orthogonal functions 
derived from a ten-level analysis of 
geopotential. 
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where the underlying tilde signifies a column vector, such as 
b19(4,9) 
9 = o9(4,8) (4.15) 
64 (4,8) 


Here A is longitude, 6 is latitude and the subscript is the 
level number. A is the Jacobian given in (4.2). 
Multiplication of (4.14) by one of the empirical orthogonal 
fume trons, e,!, gives 

mM, = Vf Vy; + A, - v 6. = 0, (4.16) 


and the variational integral imposing this constraint is 


ree + B( W,- ee + 2)m,dA (4.17) 
A 


The solution of (4.17) exactiy follows that for (4.1). 
heperes./) 1s solved for N of the most significant 
empirical orthogonal functions, the balanced fields are 


constructed from 


N 
dp = y d. E (4.18) 
~ j=] ~ 
and 
N 
7X We = 8 (VX W.)E, (4.19) 


j=] 
As shown by (4.19), the balanced vorticity is constructed 


from the empirical orthogonal functions, and then the 
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balanced winds are computed from vorticity. The degree of 
filtering is governed by the size of N, which was set to 
four to keep the jet streams sufficiently strong. 

The filtering described in this section greatly 
improved the erroneous temperature problem. A temperature 
adjustment profile for a filtered and unfiltered balance is 
shown in Fig. 12. As has been typical, the lowest layers 
were particularly poor in the nonfiltered balance, where 
adjustments were in excess of 39. Zonal averages of the 
temperature adjustments of the lowest layer in the tropics 
were -39C to -49C. The filtered balance, on the other hand, 
produced adjustments that were much smaller. The zonal 
averages of the adjustments were everywhere less than 0.59C, 
and individual adjustments in the lowest layers were about 


PC. 


Cc. WEIGHT STRUCTURE 

Typically, the weight a variable receives during the 
variational balancing depends on an intricately computed 
error structure function, which can be derived from optimum 
interpolation methods such as given by Gandin (1963). Such 
error structures are not easily derived from successive 
correction methods, but the amount of data that influences 
any point can be used to infer analysis accuracy to some 


extenit. 
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Adjustments made to the temperature in order to 
balance wind at 60°9E, 30°N on 27 Mar 1982. 

No vertical coupling (dashed line) and vertically 
coupled using empirical orthogonal functions 
(solid line). 3 
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In the methieds described in this se€tion, the weight on 
the geopotential is unity and the weight on the winds is 8, 
thereby keeping the weight function a single variable. 
Limitations to g were necessary in the region bounded by 
+309 latitude. Here, g had to be at least 107% m@ sec7* to 
avoid convergence problems. Outside this region, no such 
‘Timitations were found to exist. 

The wind weight has two basic parts. One part is 
purely a function of latitude and is prescribed without 
regard to observation density. For this part, the weight 
was set at 40,000 m¢ sec~* over the region bounded by +30° 
and at 4,000 m¢ sec~* over the remainder of the globe. This 


means that a 5m sec7l 


modification in wind would correspond 
to a l00 mor a 30 m change in height, depending on where 
the change took place. The second part of the weight 
depends on the number of observations used to describe mass 
or motion, and has arange of values from -2000 m* sec-* to 
2000 m* sec™*. For example, if a motion variable were 
influenced by five or more observations and there were no 
mass observations, the second part would have a value of 
2000 né secé, But if the reverse were true, i.e, the mass 
variables were supported by five or more observations, then 


the second part would have a value of -2000 m¢ sec7*. The 


final weight value is the sum of the first and second parts. 
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In this chapter, the classical balance equation has 
been incorporated into a global initialization that uses 
calculus of variations. The method involves the solution of 
elliptic equations on a sphere, but otherwise the method is 
simple and flexible in terms of variable weighting. This 
system was thoroughly tested by Barker (198la). In the next 
section, the role of the balance equation in data assimila- 
tion will be compared to that of the nonlinear normal mode 


method. 
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V. COMPUTATIONAL RESULTS 


Ns EXPERIMENT DESCRIPTION 

The results comparing the different initialization 
methods were obtained using the FGGE level IIa data set for 
the period between 16 Nov 1979 and 21 Nov 1979. The level 
Ila data were available on an operational basis, conse- 
quently no special processing was performed for this set. 
The observing systems available at this time are described 
by Fleming et al. (1979). 

Four different initialization methods were tested in 
data assimilation tests that covered the period of the data 
set. Two of the methods tested used some form of the 
balance equation, and the other two used the variational 
nonlinear normal mode method. 

In the two balance equation methods, the divergence was 
obtained from the forecast first guess. In the first 
method, hereafter referred to as BEl, the corrections to the 
forecast were balanced and interpolated to the model coordi- 
nate. This maintained a balance between mass and motion and 
preserved the divergence that was present in the forecast 
just prior to the time of the update. The second method, 
which shall be called BE2, balanced the total analysis 
(forecast first guess plus corrections) and then converted 


the balanced winds to nondivergent corrections. The motion 


TA 





variables were then treated as in BEl, except the mass 
variables were treated as the new values in the model. The 
main difference between the two methods is that the correc- 
tions were balanced in BEl and the updated variables were 
balanced in BE2. Both methods used the same weight assign- 
ment described in Chapter IV. 

Two variations of the normal mode approach were used. 
In the first case (NM1), the weights were varied with 
latitude. Poleward of +30 degrees, the weight on the 
geopotential was 10, that is an order of magnitude larger 
than naturally occurs in the normal mode method. Further 
increases in geopotential weight make only slight improve- 
ments in the fit of geopotential, as discussed in Chapter 
IV.0. Equatorward of +30 degrees, the geopotential weight 
was 0.5. In the second case (NM2), the geopotential weight 


was unity everywhere. 


Bs DIFFERENCES BETWEEN BALANCED AND ANALYZED VARIABLES 
The objective analysis method described above is used 
to fit the available observations relative to the first- 
guess fields from the forecast model. Unfortunately, large 
regions of the earth have inadequate data coverage in terms 
of quantity and quality. Mass corrections are frequently 
made in regions mH heUe motion corrections and vice versa. 


This, in turn, causes large imbalances between the mass and 


La. 





wind fields which often cause the model to reject much of 
the analyzed information (Daley, 1980; Daley and Puri, 
1980). Similarly, the initialization may reject updated 
information by balancing to the motion variables when the 
mass variables are more correct. To show how the balanced 
conditions differ from the analyzed ones, the differences 
between the balanced and analyzed variables were plotted for 
the different balancing methods. 

The RMS averages of the analyzed corrections and the 
differences between the balanced and unbalanced conditions 
are shown in Fig. 13 for BEl and BE2. Except for tempera- 
ture, the two methods produced about the same modifications 
moethe anaiyses. Surprisingly, these changes to the 
analyses were nearly as large as the original corrections. 

The globally averaged temperature modifications are 
much larger for the BE2 method than the BEl method, 
Primarily because of the adjustments made to the lower 
levels of the model. These adjustments were much smaller in 
subsequent update cycles ( 0.8 degrees), but were consis- 
tently twice as large as the BEl method. 

The large differences between analyzed and balanced 
winds are related to the inability of the objective analysis 
to project the corrections onto the rotational component of 


the wind. This is further illustrated in Fig. 14. 
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Regions such as the one shown have few wind observations. 
When there is an isolated report, it produces a correction 
that is largely divergent (see Barker, 1981b), and regard- 
less of the weighting on the winds, this divergence is 
rejected during the balance. The rotational component, 
however, can be fit as accurately as desired (see Chapter 
IV.D), by adjusting the wind weights. The vorticity 

Shown in Fig. 14 is reasonably close to the analyzed 

values considering the moderate values of weighting applied 
to the winds in this area, (SH)? ~ 4000 mé sec74, 

The RMS average differences between normal mode 
balanced and unbalanced variables were also compared to the 
RMS average of the analyzed corrections (see Fig. 15). In 
these comparisons, however, the divergent component of the 
wind correction was removed from the analysis before inter- 
polation to the model coordinates. The method with stronger 
weight on geopotential produced closer fits of temperature 
and pressure than did the other method, but it failed to fit 
the winds as well. This is to be expected in the varia- 
tional method. The globally averaged wind differences were 
between 1.0 and 1.5 m secqt, but when divergence was not 
removed from the analysis, the differences were consistently 
larger (between 2.0 and 2.7 sec7t), The surface pressure 


differences were as large as the analyzed corrections over 


much of the earth. 


8 | 





Comparing the balance equation results (Fig. 13) to 
those of normal mode methods (Fig. 15) indicates that both 
techniques were similar in the way that they modified 
analyzed variables. The balanced winds differ from the 
analyzed winds by about the same amount for all methods when 
the analyzed divergence was not removed prior to balancing 
(these curves are not shown). The normal mode methods were 
slightly better at fitting the analyzed temperature, but 
poorer at fitting the surface pressure. 

The modifications to the analyzed variables were 
Significantly large in all experiments performed. However, 
whether the impact of the balancing can be considered 
beneficial to the data assimilation process principally 
depends on the magnitude of the gravity wave noise that 
still exists. In the next section, the noise levels that 
are present in the model before and after balancing are 


compared for the balance equation and normal mode methods. 


Ge ELIMINATION OF GRAVITY WAVE NOISE 

In this section, gravity noise produced in forecasts 
initialized with no balance, the balance equation and the 
normal mode method are compared. In these comparisons, 
surface pressure tendency was used as a measure of the 
gravity noise. Although this quantity reflects only the 


presence of external gravity waves, the highest frequencies 
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come from these external waves. Except where stated 
otherwise, adry version of the global model was used to 
perform the integrations. 

In the first comparison, the effectiveness of the 
balance equation was tested by comparing it with an 
initialization that just removed divergent winds from the 
update corrections. In the balance equation method used for 
testing, the corrections were balanced, interpolated to 
model coordinates and then added to the forecast first guess 
(method BE2). As can be seen in Fig. 16, the balance equa- 
tion forecast is less noisy than the forecast initialized 
with nondivergent winds, particularly during the final hours 
of the forecast. The slower adjustment rate of the 
unbalanced forecast is probably due to the scale of the 
gravity waves, since global models may carry large gravity 
waves that do not readily disperse their energy (Bourke, 
1972). This version of the balance equation initialization, 
on the other hand, is effective at removing these large 
scale waves, but it does not balance around terrain. 
Consequently, the balance equation method contained notice- 
able small scale noise. | 

Using the model's normal modes, the linear balance was 
compared to the nonlinear balance (Fig. 17). The nonlinear 
balance forecast required one hour adjustment time and 


resulted in about one half as much noise as did the linear 
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balance forecast. These results illustrate the effective- 
ness of the nonlinear balance step, which can balance around 
terrain (Daley, 1979; Williamson and Temperton, 1981) and 
can include the nonlinear component of the balance 
(Machenhauer, 1977). 

Plots of surface pressure tendencies from predictions 
initialized with the balance equation and with the nonlinear 
normal mode method are shown in Fig. 18. These results show 


the superiority of the normal mode method over the balance 


equation. After seven hours, however, the forecast started 


from the balance equation contained only slightly more noise 


than the forecast from the normal mode method. 

Finally, the impact of adding latent heating effects to 
the forecasts is illustrated in Figs. 19 through 21. fhe 
increased noise level is most noticeable in the normal mode 
runs (Fig. 20) where the forecast seems to have required 
about six hours of adjustment time. The balance equation 
run was only slightly affected by the heating (Fig. 19), 
possibly because less precipitation was generated during the 
early hours of the forecast. In any case, the result is 
that the normal mode method was no better at noise 
Suppression than the balance equation when jiatent heating 
effects were included in the forecast. This result was 
observed consistently several times during the course of 


development of the methods. 


ad) 





";apow yseddVOJ BUY UL POpN_DUL You st BHutzeay 


yuayey, uaym (aut, pLtos) votyenbsa aoueyteq ayy wos pue (auUL| paysep) 


UOLJEZLFeLALUL 


at 





apow j;ewuou ayy wo} ADuapuay aunssaud adesuns JO SWYH ‘BT 


YH NI NV 
oT 8 9 a © a 

vea* 
val BG) s 
=e 
W 
THR’? + 
Sen = 
m 
Ww 
Ww 
= 
= 
foe’ y 
-ZAQ’ _, 
= 
Ze 
o 
m 
z 
o 
~< 

caa*~ 
1 : 
ow 
™ 
W 
m 
O 

FAR’ 


u) 
LOR 
' ale 


aayth 


88 








“out{ ptytos e Aq paquasaudau 


St $299sJa JeOY Yuazey, BHButpnypout ysedauos |JYY pue aUL{ paysep e 


Aq paquasdudau St 
UVLM Pazt{etylut 





ySeIIVOJ DLYeQetpe auy “pouzyow uotyenba aogue{eq ayy 
Sysedau0s ut AJQuapusay aunssaud adsejsuns jo abeuane swy 


YH NE AVL 





“6T 


a6B8° 


‘OD 

‘Se 

le 
® 


iB a) 
bot 
he) 

6 


lanl 


“RONSGNSL 3YNSS3uad SWY 


33S /aW 


89 








-anbluyoay apow ;ewsou seauL~UOU 3yy AOJ 3d]99KI GT bite O 7 ele i, tans 


YH NI AVL 





Siz °6 


600° 


788° 


a 


S6a° 


L4 


AINSONSL SYNSSAYd SWY 


¢ 


J4S/9W 


90 





"3S5e994U0J JU UL papnNy_ Dut St Hulzyeay yuaze,_ UDYM 
(QuL{ pl_os) uolyenba aduejyeq ayy soy pue (auUl| paysep) uolzezipelyztul 


apow jewu4ou ueaulL{uoU ayy sos ADUSpUdy JuNnssadd adejsuns JO SWY 
YH NI NVI 
rz | of 9 y A iy 
eag° 
¢Ga° 
PQ” 
269° 


6Z26T AON 9T LWD 00 





— 86a" 


"Iz ‘bly 


“AONSQNSL 3SYNSS3ad SWY 


93S /dW 


S) 





In summary, the normal mode method was most affected by 
latent heating in the forecast. This was probably because 
it allowed more precipitation to occur during the early 
hours of the forecast compared to the balance equation 
method. Considering that one problem associated with 
initialization has been the lack of precipitation early in 
the forecast, the noise level increase in the forecast 
including latent heating may be a symptom of a beneficial 
result. It is noteworthy that latent heating effects 
generated twice as much gravity noise as the integrations 
that did not include latent heating. This was true even 
after the initial adjustment period was complete. To 
explore the question of latent heating effects more fully, 
comparisons of the precipitation, vertical motion and 
cyclogenesis during data assimilation runs were made. The 
results from these runs are described in the following 


section. 


Dy. VERTICAL MOTION PRECIPITATION AND CYCLOGENESIS 

The lack of precipitation early in the forecast is 
considered to be a major problem in the initialization of 
numerical models (Tarbell el al., 1981; Bengsston, 1981; 
waylevaeel9el). this is particularly true of forecasts 
initialized without vertical motion. An example of the 


problem is shown in Fig. 22. These two forecasts were 
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initialized with the balance equation with no provisions to 
paciuge wertical motion. The forecasts for the region 
west-northwest of the Hawaiian Islands predicted large 
precipitation rates after six hours, but only slight amounts 
before this time. 

These initially small precipitation rates occur when no 
vertical motion is included in the initialized data. Inclu- 
sion of the omega equation (Tarbell et al., 1981), nonlinear 
normal mode initialization (Leith, 1980) and retention of 
the forecast first-guess vertical motion provide possible 
solutions to this problem. Unfortunately, the omega equa- 
tion and normal-mode methods have not worked well in the 
erepgacs (BeWgsston, 1981; Tribbia, 1981), and the forecast 
first guess may suffer from inaccuracies in the forecast. 

The balance equation methods discussed in Chapter [II 
use the forecast first guess divergence. The normal mode 
initialization partially recomputes this divergence through 
the nonlinear balancing of the external and first internal 
vertical modes. The forecast and derived divergence from 
the normal mode methods are given in Fig. 23. A sequence of 
forecast and computed divergence fields is given in Fig. 24. 
Notice that only small differences exist between the fore- 
cast and computed divergence patterns. inis similarity 
between the two divergent winds suggests that the forecast 


Sivermgence’is a fairly accurate quantity. 
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Forecast first-guess divergence (a, c, e) and 
divergence computed using the normal mode method, 
NM2 (b, d, f), for three successive 1l2-hourly 
updates beginning at 12 GMT 1/7 Nov 1979. Contour 
interval is 1072 sect. _The interval between 
-1-1072 sect and -2-°107? sec7l is cross-hatched 
and the interval between 1075 sec~+ and 2-10~° 
sec™+ is cross-hatched twice. 
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Precipitation rates from five forecasts that are 
identical except for the initialization procedure are shown 
miamge 25. Im the first case (Fig. 25a), no initialization 
was used. Even the excessive divergence produced by the 
analysis was included in the initial data. The other 
methods, which are designated BE1, BE2, NM1 and NM2, are 
described in Chapter V.A. These forecasts all produced a 
nearly identical precipitation pattern along the coast of 
North America. The normal mode methods (NM1, NM2), however, 
produced less precipitation in the central Pacific than did 
the balance equation methods (BE1, BE2). However, on the 
Other hand, normal mode methods have produced the largest 
precipitation rates south of Japan. None of the methods 
produced a persistent bias in the strength of precipitation 
paees. Considering the sensitivity of precipitation rate to 
small changes in initial data, however, the similarity of 
the patterns is quite remarkable. In particular, the lack 
ef spurious precipitation in the forecast without any 
balancing is most surprising. 

Tie: precipicatien rates durigag the first twelve hours 
of forecasts initialized by the balance equation and normal 
mode methods are snown in Fig. 26 for a region just north of 
South America. Separate rates are given for the first and 
second six-hour periods to illustrate the impact of initial- 


fame VOurdng the first six hours of the forecast, the 
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precipitation extended over much larger regions and was 
stronger than during the second six hours. This tendency is 
just the opposite effect observed in the nondivergent case, 
and helps to explain why the model produces more gravity 
wave noise when the latent heating effects are included in 
the forecast. The precipitation patterns appear to be less 
Similar between the forecasts during the second six hours 
tian during the first six. In fact, after several update 
cycles, the precipitation fields produced by the different 
methods contain little similarity in the tropics. 

Slight differences due to the various initialization 
methods may cause the assimilation runs to diverge slowly 
from each other. The largest differences are likely to 
occur in regions of strong baroclinic instability where 
precipitation cam play a role. To examine this possibility, 
a moderately intense surface low, which developed along the 
Aleutian Islands, was studied. This case of surface cyclo- 
genesis of 20 mb in 24 hours was supported by an upper level 
Short wave that traveled along the island chain in twenty- 
four hours. Three twelve-hour forecasts are shown for the 
various initialization methods tested (Figs. 27 through 34). 
Comparing forecasts rather than analyses nelps to eliminate 
the possibility that gravity modes allow a closer fit to the 
data than actually exists by the meteorological modes. As 


Can be seen from the verifying analyses (Fig. 35), none of 
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Twelve hour forecasts of 500 mb geopotential (a, 
c, e) and sea-level pressure (b, d, f) during the 
assimilation run using the BEl initialization 
method. Contour interval for geopotential is 60 m 
and for sea-level pressure is 4 mb. The 4920 to 
4980 m interval is cross-hatched in the 500 mb 
maps. The 996 to 1000 mb and 980-984 mb 

intervals are cross-hatched in the sea-level 
pressure maps. 
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the forecasts developed the system as rapidly as the 
atmosphere. However, the normal mode method came closest by 
developing slightly deeper systems each time. The NM2 
method managed to do slightly better than the NM1 method. 
The weakest system was produced by the BE2 method, but it 
was only 4 mb too weak. The normal mode methods generated 
the largest precipitation amounts. This is particularly 
Erwme dur wing thre 6 to 12 hour forecast*made from the 17 
November 1979 12 GMT data. During most of the period, 
however, the precipitation differences were quite small. 

In general, forecast experiments that used divergence 
in the initial conditions produced precipitation rates in 
the early hours of the forecast that were larger than those 
produced after six hours. This is the opposite effect 
observed in forecasts using a nondivergent initialization. 
Whether the divergence was computed from the normal mode 
method or derived from the forecast first-gquess seemed to 
make little difference in the resulting precipitation and 
divergence. During data assimilation experiments lasting 
several days, however, these differences detween the methods 
became more pronounced. A study of cyclogenesis over the 
North Pacific showed that while the normal mode methods 
tended to deepen the system faster and somewhat more 
accurately, the balance equation methods were nearly as 


effective. 
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Unfortunately, these results which describe a single 
case are hardly conclusive. To extend the verification of 
the different initialization methods over broader areas and 
more cases, comparisons of many forecasts to observations 


are made in the next section. 


on FORECAST VERIFICATION 

Justification for balancing during data assimilation 
runs is readily demonstrated with the verification compari- 
son of a balanced and unbalanced forecast. Forecast 
verifications against all observations for these runs are 
plotted versus latitude in Fig. 36 for geopotential and 
winds. The results indicate that although the wind 
verification is unaffected by the presence of gravity wave 
noise, the geopotential verification is drastically 
oecected. Mis forecast from an unbalanced initial state 
has RMS errors almost double that of the balanced forecast 
over much of the earth. Such large errors in a forecast 
first guess may cause the quality control programs to reject 
observations that should not be rejected and to add noise to 
the resulting analysis. (Heavy damping filters applied 
during integration may make the forecasts more usable, 
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Latitudinal variation of RMS (a) pressure height 

and (b) wind differences between observations and 
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While it is certain that the forecasts from an 
unbalanced initial state are unsuitable as a first guess for 
an analysis of geopotential observations, and therefore are 
not useful in data assimilation, the most effective 
balancing procedure to be used is not clear. Some small 
differences do exist between the assimilation results using 
different initialization methods. This can be seen (Fig. 
37) in the analysis of a short wave upstream of Australia. 
The 12-hour forecasts (rather than analyses) are shown for 
this system so that data resolved onto gravity waves could 
be dispersed by the model. The normal mode methods (NM1 and 
NM2) produced a slightly deeper wave than did the balance 
equation methods (BEl and BE2). This is consistent with the 
nesmlcs in Chapter V.C. 

An extensive objective verification study was performed 
for the different methods, where the RMS differences between 
the 12-hour forecasts and observations were computed. These 
computations were made for ten forecasts covering the period 
of the data assimilation experiments. Unfortunately, this 
type of test may tend to favor the smoother forecasts, and 
therefore the interpretation of results should be made with 
this in mind. 

Although verification was performed by region as well 
as globally, the regional verification added no new informa- 


tion not readily apparent in the global statistics. The 
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verification is presented by observation type, which tends 
to regionalize the verification to some extent. For 
example, observations from ships, satellites and aircraft 
are primarily taken over the oceans, whereas radiosonde 
observations are mostly taken over land. 

Forecast verification against radiosonde geopotential 
and wind observations are presented in Figs. 38 and 39 for 
the four methods of initialization. The lowest RMS differ- 
ences between forecasts and observations occurred in the 
balance equation method, BE2. However, the differences 
between methods are very small. For geopotential, the 
differences produced by the various methods are generally 
less than 2m, which is only about 4% of the total RMS error. 
For wind, this difference is generally less than 3 m sec7l, 
which also represents about 4% of the total error. 

Comparing the normal mode methods, it can be seen that NM1 
verified against geopotential observations slightly better, 
whereas NM2 verified against wind observations slightly 
better. This reflects the emphasis that the variational 
Dalancing placed on the respective variables. However, 
notice that BEl verified geopotential as well as the NM1 
method and wind as well as the NM2 method. 

Forecast verifications against ship observations of sea 
level pressure, which were converted to geopotential through 


the hydrostatic equation, are given in Fig. 40. The 
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improvement of the BE2 method over the other methods is much 
more significant for this data. Compared to BE1, the 
improvement is as high as 10% of the total RMS error. The 
NM1 method, which weighted the geopotential analysis more 
than the NM2 method, verified somewhere between BE2 and BEl, 
whereas NM2 produced approximately the same results as BEl. 

Verification using satellite-generated geopotentials 
(Fig. 41) again show the BE2 method to be superior to the 
Other methods. The improvement ranges between 6 and 12%. 
The other initialization methods resulted in rather similar 
verifications. 

Satellite wind forecast verification (Fig. 42) shows 
that the BE2 method gave about 0.5 m sec} smaller forecast 
error, or about a 5% improvement over the other methods. 
This verification is based mainly on the low level ( 925 mb) 
satellite observations in the tropics. Once again, the 
Other verifications are similar to each other. 

Unlike the other types of data, aircraft wind forecast 
verification failed to show much difference between the 
methods (see Fig. 43). Since aircraft observations are 
primarily taken between 300 and 250 mb, this suggests that 
the four methods produced comparable results at the upper 


levels over the oceans. 
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In summary, the BE2 method, which uses the balance 
equation to variationally balance the complete analysis 
(forecast first guess plus corrections), produced the best 
verification scores. This result is most pronounced at the 
low levels over the oceans. Visual inspection of the fore- 
casts produced by this method indicate that they were also 
smoother than forecasts produced by other methods. While 
RMS scores tend to be better for smooth fields, the much 
better verifications against ship observations produced by 
the BE2 method are difficult to dispute. The differences 
among the remaining methods are too small to conclude which 


is best. 
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VI. SUMMARY AND CONCLUSIONS 


In this study, two different static initialization 
methods have been developed and tested in data assimilation 
experiments. These methods not only reduce gravity wave 
noise, but they also fit the meteorological modes to the 
most accurately analyzed variables through calculus of 
variations. The methods are based on the balance equation 
and a nonlinear normal mode procedure. Of the two methods, 
the normal mode method would be the more cumbersome to 
apply, because the most general form with the weighting 
fully variable in the horizontal requires more computer 
memory than was available. 

The two methods were compared in various ways. The 
comparison tests were designed to show how the analyzed 
fields were modified during the balancing, how each 
controlled gravity wave noise, how precipitation varied 
during the early hours of the forecast, how the forecast of 
a rapidly developing system was affected, and finally, how 
short range forecasts from the various methods verified 
against observations. 

Both methods modified the analyzed fields during the 
baiance by a large amount when compared to the size of the 
corrections (differences between analysis and forecast first 


guess). Surface pressure modifications were particularly 


176 





large, sometimes exceeding the magnitude of the corrections. 
Modifications to the winds were primarily caused by 
unrealistic divergence patterns in the analyzed winds, which 
were nearly as large as the analyzed corrections. When the 
divergence was removed prior to balancing, the modifications 
were much smaller. For wind and surface pressure, the 
magnitude of the adjustments was very nearly the same for 
the two methods. However, the modifications to the tempera- 
ture analysis were larger for the balance equation method 
when it was used to balance the total analysis (corrections 
plus first guess) than when it was used to balance just the 
Corrections. 

In terms of gravity wave noise removal, the normal mode 
methods performed significantly better in a dry version of 
the prediction model than did the balance equation method. 
Adding the effects of latent heating, however, tended to 
mask this improvement, because the heating effects caused 
gravity wave noise to more than double relative to the dry 
version of the model. 

The increased gravity noise due to tne heating may be 
caused Dy many factors. It is partly caused dy the way in 
which the heating effects are included. Each fifth time 
Step, the heating effects are added. A Matsuno time step is 
then used prior to resuming the leapfrog time stepping. The 


gravity wave amplitudes produced are larger than if the 
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heating effects are added incrementally over each time step. 
Another factor that may make the latent heating effects more 
pronounced is that the model tends to be warmer and drier in 
the lower layers than observed in the atmosphere (Johnson, 
1976; Payne, 1981). After each update, the solution tends 
toward the model equilibrium state. Also, because no 
analysis is made of moisture, changes in the mass and motion 
structure may require that the model undergo some adjustment 
before latent heating matches the updated systems. None of 
these effects are controllable by the initialization proce- 
dures under study, however. 

Integrations performed without divergence (vertical 
motion) required several hours to develop realistic 
precipitation rates. This was particularly true in the 
tropical regions in a test with the balance equation 
initialization that did not include divergence from the 
forecast first-guess. 

Comparisons of divergence from the forecast first-guess 
and that computed using the nonlinear normal mode balance 
applied to the external and first internal vertical modes 
revealed only small differences. Other comparisons of 
precipitation forecasts from the balance equation, normal 
mode and no balance conditions showed oniy minor differ- 
ences. In the balance equation and no balance methods, the 


forecast first-guess divergence was present in the initial 
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conditions. The unbalanced forecast surprisingly produced 
no noticeable spurious precipitation forecasts, even when 
the unrealistic divergence produced during the wind analysis 
was included. Unlike the forecast initialized without 
divergence, however, largest precipitation rates were pre- 
dicted during the first few hours of the forecast. It 
appears, then, that whether the forecast first-guess 
divergence is partially recomputed using the normal mode 
methods or even mixed in with unrealistic divergence, little 
difference will exist in the precipitation forecast. 

The effectiveness of the model in assimilating data 
around a rapidly developing surface low over the Pacific was 
examined for the different initialization methods. The 
representation of the cyclone development was very similar 
for the various methods. The maximum central pressure 
difference between the methods was 4 mb. The more intense, 
and slightly more accurate representation of the low was 
produced by a normal mode method, whereas the least accurate 
was produced by the balance equation method that balanced 
the total analysis and used the forecast first-guess diver- 
gence. The results from this single example are only 
suggestive, however. 

To produce results over a wider number of cases, the 
forecast first-guess fields were verified for the four data 


assimilation methods. Two methods used variations of the 
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balance equation method and two others used variations of 
the normal mode method. The results from these runs showed 
that while minor differences existed between the forecasts 
for much of the data, the balance equation method that 
balanced the total analysis, rather than the corrections, 
produced the most accurate short-range forecasts near the 
Surface over the oceans. The forecasts produced from the 
normal mode methods tended to contain slightly stronger 
systems than the balance equation methods, which may have 
actually reduced all the verification scores for this method 
relative to the smoother fields from the balance equation 
method. 

In conclusion, the results of this study show that it 
is possible to initialize a model with the forecast first- 
guess divergence. This allows continuity in the precipita- 
tion rates produced by the model during the updates. 
Consequently the variational balance equation method 
produced results that are competitive with, and in some 
respects, better than the more elegant nonlinear normal mode 
method. This conclusion is based on precipitation rates 
during the early hours of the forecast, gravity wave noise 
and short-term forecast results. In terms of variational 
weight assignment, the balance equation methods are less 


cumbersome, and thereby allow more flexibility. 
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As a note of caution, however, it should be mentioned 
that assimilation systems, such as the ones tested, are too 
complicated to guarantee absolutely that the tests were 
without flaws. For example, the balance equation method 
that scored the highest forecast verification scores was 
also the only method that used a slight variation of the 
interpolation to model coordinates. It interpolated 
updated mass variables instead of corrections. Interpola- 
tion of corrections did not insure that the sea level 
pressure under terrain matched observations corrected to sea 
level. Consequently, regions such as the Himalayas 
contained sea-level features that were not present at 
terrain level. Another factor is the type of assimilation- 
prediction that was used. The version of the UCLA model 
used in the assimilation produces much larger surface 
pressure tendencies than does a dry version of the same 
model. This factor, which masked the benefits of the normal 
mode initialization, may not be so prevalent in future 


generation models. 
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APPENDIX A 
DATA ANALYSIS WITH SIMULTANEOUS FILTERING 


The inherent smoothing of a successive correction 
scheme can be determined by neglecting discrete spacing of 
the observations and assuming isotropy in the weight 
Boeeificanion (barmes, 1973). Under these conditions, 

F(r) -{*, (e)f(r-e)de. (Al) 
The convolution theorem allows us to take the Fourier 
transform of this equation or 
F(K) = 3(K)F(K) (A2) 
where K is the wave number equal to 27/i, and A’ iS wave- 
length. The hat is used to show that the variable is 
transformed into wave number space, e.g., 

ue) ee eee ode (A3) 


Recall that F is a complex quantity whose magnitude is 


represented by 


i} 


|F| = VFF* (A4) 


or 

Le | 3c Rea (A5) 
which means that the spectral response of the analysis may 
be identified by |a(K)}. 


If the analysis is repeated to converge on the data 


more closely, the resulting product will be 


tod 





soclule  I2 we f° wate )[f(r-c)-Fy(r) ]de (A6) 


where Fy is the result of applying (Al) with weight function 
wy. This equation transforms to 
Fr(K) = Fy(K) + do(K)F(K) = do(K)FY(K), (A7) 
or 

Fr(K) = [az(K) + d9(K) - dp(K) ay (K)IF(K) (A8) 
and the spectral response of two passes is identified by 

|w | = | wy : wo WoW | (AQ) 


Specifying the weight function as Gaussian, where 


wy (e) = a) exp(-2°/B°) , (A10) 


i 


wole) = ay exp(-e°/B°) , (Al1) 


Simplifies the computation of spectral response. y and Be 
are arbitrary constants and a, and a> are normalizing 
Coenrricivents, e.g., 

a, = if exp(-e2/8¢)de]7t. (Al2) 
The Fourier transform of (Al0) and (All) are 


“aA 


wy (K) exp (-B°K°/4) and CAalsn) 


exp(-yB°Ke/4) . | (A14) 


i} 


w(K) 


It is now possible to design the spectral response, wy, dy 


appropriate choices for y and 8B. 
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The objective analysis for discrete spacing of the 


observations is 


Pane) oes) t ineiel) ey) 


bo 
al 


hietums equation and those tnat follow in this section, the 
upper case letters, such as F', designate gridded fields of 
the variable being analyzed, whereas the lower case letters 
are related to observations that are irregularly located in 
space and time. The vector, e,, specifies the separation in 
space and time between the observation, fiend the ‘gidieat 
point r. The quantity f' represents the difference between 
the observation and the forecast that is to be updated and 
mois melemaindiyzedscomnection field that results from Al5. 

The weighting function, y, determines how the observa- 
tions are weight-averaged at each point on the grid. There 
mumo Wooer jimit to tie size of «, but a practical limit 
CeeuinS when ay 1S sufficiently small as to make correction 
meaningless. fhe volume over which computation is performed 
is referred to as the scan volume, and it represents the 
region from which observations are weight averaged to 
produce a value for point r. eauataon AlS@was first 
developed by Bergthorsson and Doos (1955). Cressman (1959) 
overcame some of the inherent oversmoothing of this 


technique by incorporating multiple analysis scans, with 





each scan using a progressively smaller radius. [In this 
Successive correction metnod, the second pass uses a weight, 


w2, to reduce the remaining discrepancy between the analysis 


after the first pass and the observation, f", i.e., 


F"(r) = wole,p)f"(r-ey). (A16) 


Hom 


k=] 


ines resulting GCorrection is equal to the sum of F' and F". 
Although the Cressman method is usually designed to 

perform three or more passes through the data, Barnes (1964, 
1973) and Stephens (1967) have shown that careful design of 
the weight function may provide adequate results after only 
two passes. The weight functionsdefined by Barnes are 

Wi 7 a4 exp(-¢ ¢/B¢) (A17) 
and 

MQ = ey, exp(-« ¢/yB*). (A18) 
As already shown, an analysis that uses (A118) and 
(Al19) to compute the weights, and has a sufficiently dense 
observational coverage, will produce a spectral response 
Or(K) = exp(-B¢K*/4) + exp( - BCK*/4) 

exp[- (1+y)8°K°/4], (A19) 

where 3 and , are arbitrary parameters, and K is the wave 
number (2,/,, where \ is wavelength). The symbol (*) 
Signifies that the variable is transformed into wave-number 


Space. It is possible to use this equation to design 





variable filtering characteristics that are dependent on 
observation type in the analysis. The parameter B is used 
to limit the volume of influence of an observation, and the 
parameter , is used to specify the degree of inherent 
famecer ing . 

To design the desired filtering into the analysis, it 
is convenient to define the weight function as the product 
Or ccunmee functions, i.e., 

w = upl(en)woy(ey)wzelez) (A20) 
where the subscripts H, v and t designate the functions 
describing the horizontal, vertical and temporal dependence, 
respectively. In this form, it is possible to make the 
filtering different for each dimension in the analysis. 

A vertical weighting function that allows ample 
vertical variability while maintaining some vertical 
coupling is desirable. The weight function, Dy. ine AZO) 
is proportional to exp(-c,/B,°), Wine Mge cr meen mt) Pie.) 
represents the pressure separation between the observation 
and analysis level. The constant 8, is 0.6, which produces 
a vertical scan radius that corresponds to the positive 
values of the prediction error covariances computed by 
Hallett (see Rutherford, 1976). The vertical filtering 
parameter, y,, is allowed to vary with observation type. 
For satellite soundings it is 0.8, whereas for other data 


types it is 0.3. Thus the satellite sounding corrections 





are very smooth with height, thereby avoiding the problem of 
removing inversions that are formed by the model or are 
present in the other observations (Tracton et al., 1980). 
For example, updates from satellite soundings containing 
vertical wavelengths of 0.5 are damped to 20% of their 
original values, whereas this same wavelength in a radio- 
sonde is damped only 50%. 

As is the case for the vertical function, the 
horizontal weighting function is designed to limit the 
influence of an observation to an area that roughly 
corresponds to the positive values of correlation structure 
function. Using the curves produced by Buell (1972), By is 
then 3.24 grid intervals on the 2.59 mesh used by the 
analysis. To account for the spherical geometry, the hori- 
zontal separation between observation and grid point is 
defined by 

eye = (ux)? + y?, A2ay) 
where x and y are the zonal and meridional distances in grid 
intervals, uis the map factor 

f= max [eoss, 0.95], (A22) 
and 9 is latitude. The lower limit on y distorts the region 
of influence an observation may nave poleward of 609, but it 


prevents obvious computational probiems near the poles. 





Inman (1970) has shown that successive correction 
methods have a tendency to diffuse narrow jets such as occur 
in the upper troposphere. An ellipsoidally shaped weighting 
function with the major axis aligned along the wind 
direction tends to avoid the difficulty. Therefore, the 


weighting function for the wind is computed from 


wiley) = exp(-e4 m./Bu°) » (A23) 
where 

ye = On] +. 3.40 SRK Cea) Hees (A24) 

7 arctan (v/u) and (NZ oe) 

9 = arctan (y/Xx). (A26) 


Inclusion of the map factor in (A26) made no impact inthe 
analysis and therefore was omitted to save considerable 
computer time. 

The present assimilation system updates the forecast 
every twelve hours, which means that during any single 
analysis time, there are likely to be data whose observation 
Cime differs from that of the analysis by six hours. As is 
the case in spatial dimensions, a poorly defined time weight 
function will result in damping and shifting of the smal! 
scale waves. If the weight in time is determined from (Al7) 
and (A183) where <;} is the separation in hours between 
observation time and analysis time, the inherent tempora|] 
filtering of the analysis depends on the values of 8B; and 


Y + Pewiicedcisdwamme to compute 8+ from time correlations of 


ies 





different data types, however very little literature exists 
on this type of study. Barnes (1973) used the phase speed of 
the meteorological systems to compute values for B; that are 
consistent with the spatial parameters. Assuming a phase 
speed, C, of the order of 20 kt, then By can be determined 
from By/C, which is about 20 hours. For the update interval 
used, however, the temporal variations of the weights is 

only about 10%. 

The amount of horizontal and temporal filtering 
infverent in the analysis is governed by yy and y;>, 
respectively. The observational accuracy may be a factor in 
the determination of these values, as observations 
containing large random error tend to produce fictitious 
Short waves which should be filtered more than the longer 
waves. The density of reports also becomes a factor, since 
aliasing may result when insufficient reports are used to 
describe a wave (Stephens, 1972). The filtering parameters 
selected (YH ¢ = 0.3) give a response for the four-grid 
increment wave of only 25% with a fairly steep rise to 80% 
for the eight-grid increment waves. 

The error checking procedures used in analysis schemes 
may account for sizeable differences bdetween various 
techniques. Even the most sophisticated schemes are not 
immune to large errors in the analysis that are caused oy 


improper rejection of the observations. The most difficult 
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problem is avoiding the rejection of observations that 
result from a poor forecast. This is most likely to happen 
if the forecast is used to check the data, rather than the 
more desirable “buddy check" method, in which observations 
are compared. In the design of the error rejection 
procedure, an attempt was made to retain as many 
observations as possible, even at the expense of accepting 
data with errors. These erroneous data cause small scale 
effects that tend to diminish during the balancing 
procedure. 

Three separate procedures are used to detect erroneous 
data. First, the radiosonde data are subjected to a 
vertical consistency check that requires the lapse rate be 
stable. Data is corrected through interpolation from 
adjacent levels with good data whenever possible. Secondly, 
gross errors are removed by rejecting the data that disagree 
with the forecast by more than five standard deviations of 
the expected error at that level. Finally, the remaining 
observations are used to perform a single pass, two- 
dimensional analysis. Each observation is then checked by 
first removing its effect from the analysis. This is done 
by determining the effect the observation in question has on 


the analysis at the closest grid point, 
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where ue is the weight given fp at this location. The 
analysis at this point excluding the observation being 
tested becomes 

Fe Fe ope (A28) 
The observation is rejected whenever disagreement between F. 
amid) fipmexceeds three standard deviations of the forecast 
error and there exists the equivalent of two other observa- 
tions within one grid interval of the observation in 
question. This procedure rejects less than 0.2% of the 
total observations. An advantage of this method is that it 
utilizes computer code that is used to perform the analysis, 
which decreases the size and complexity of the computer 
program. 

The successive corrections procedure described above is 
both simple and very fast compared to other methods, such as 
the multivariate optimum interpolation method (Schlatter, 
1975). Furthermore, in experiments performed by Otto- 
3liesner et al. (1977), the more sophisticated and time 
consuming methods did not produce significantly better 
results. Unfortunately, the successive corrections method 
is univariate, which means that no attempt is made to 
constrain the corrections to be consistent dynamically. 

This makes the balancing component of the assimilation 
System critical if the analyses are to be optimally combined 


to produce meteorologically consistent updates. 
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APPENDIX B 
LINEARIZED HYDROSTATIC, THERMODYNAMIC 
AND CONTINUITY EQUATIONS 
The vertical grid (see Fig. Bl) has all variables except 
Pressure are defined at odd Jevels and therefore 


~ 


Pk - Pk +2 = C Pa ome ya For k = | AGS earn oo (B1) 








p 
and 
ce = es CAP Py Oy (B2) 
Interpolation formulas 
“| age “1+K 
7 a es rae 
| +c ] +k 
oe ] 1 Ps ; Peay (B4) 
k oS [+k p = p z 
0 S k=] 
and 
241 = Agape F Bear okee (Be) 


are used to produce energetically consistent equations where 


Py, = Sy eaves (B6) 


oe e (37) 
K ae 
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Aj z Paty 7 PK 
om Mean ~ Py 
. for k=1,3...k-4 (B38) 
3g ke Pa 
RT Pee 7 Pe 
awe . 
= ae ee er. 
Py 
AK.) 2 SCT. (B9) 
Pap 
K K-2 
MaKe? 
as Po 
B = - (B10) 
X=] ee: | 
: K K-72 


T is temperature, p 


@erconstant pressure. 


Peg e | 
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is pressure and C, 


is the specific heat 
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Grid configuration. 





The geopotential at each level is computed through 


integration, or 


by =o, + C, (PS - PY), (B11) 
and 
DS Og mi oon Brae.) a 
(B12) 
K 
aoe Cyt Py ~ ea) Weary, 


The primed sum indicates increments of 2. Now we can write 


the hydrostatic equation as 


P Ky 
n=] 
where 
0 n<k 
CAP ne = POAne ey 86 nek 
ee = 
kn 
CAP ee : Py inn Ula : 


> - ¢ a : (B14) 


The finite difference form of the thermodynamic 
equation (Eq. 299 in Arakawa and Lamb 1977) in orthogonal 


curvilinear coordinates is 
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» F = rua », Gs nye= , S = Io, the overbar 
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where Il = an 


1s a linear average in the direction of the variable 
indicated, and Sy is a difference taken in the sGinrece1e no 1 
the subscript. 
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6 e 
k oS) |e 
Peeeeleoe § 6G & —— i = 0, (B16) 
which gives 
~ 4K 
BP x auc eee 
k | pk K g a) 2 
1S, VN, + A P;;6,t8 bes = Negras (B17) 
1 k on _ 
om (noo) ss o¢ > 97); 
or 
7 5 T 5 
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where Ti is the rest-state temperature. 
Substituting the linearized form of the continuity 


equation (Eqs. 166-167 in Arakawa and Lamb 1972) give 


K K 
oe = ‘, (9 Ve) Ao. 1 Opes ee 7° ( V, a9, Oh 
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k=] K 
ie mee iveV) Ag. + Opi, (9 V) 45 + Q. (B20) 
n=] n=] 
6,1nn = - Bes V) Aa. + a (B21) 
and 
Vy, dP, 
(naa), = C, PL (ia? (B22) 
gives 
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The formula 


(B26) 
echo Veneer) + By ola, _,(P)-P, 5) 4 
T = = DIALECT REEEEEa tc k<K 
Ady 
Pp -p for k=K 


is that derived by Arakawa from the interpolation formulas 
(B3)-(B10). 
Note that the continuity equation may also be written in 


matrix form by defining 


AG, 
Ad>5 
t= . : (B27) 
Ady 
® 
so that 
‘. line ee (B28) 
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APPENDIX C 
EMPIRICAL ORTHOGONAL FUNCTIONS 


Holstrom (1963) showed that geopotential, o(x,y,p), 
may be accurately described in terms of a series in ortho- 


normal basis functions, o,(p), derived so that a partial 


sum, 
n ~ 
Dent Say aie) 3 “24 ag(XsY)oK(P) CEN | 
where 
Po 
ar (x,y) o(X,y¥,P)o,(p)dp, (C2) 
] 


produces an optimum fit for all choices of n. In these 
equations, the horizontal coordinates are x and y and the 
Vomredie@cvordindtee 1S p. Vie valWwes p; and p> are the 
vertical boundaries of the domain and the areal mean value 
of » at each level has been removed. 

Obukhov's (1960) method, which is computationally 
easier to use than Holstrom's (1963), uses the auto- 
covariance as a characteristic measure for determining the 
empirical orthogonal function, iji.e., 

B(p.p') = o(x,y,p}o(x,y.p'), (C3) 
where the overbar operator designates a horizontal average. 
This function describes the variance of » when p is equal to 


p', and it describes the covariance otherwise. The 








redundancy of the o-profiles is identifiable by the size of 
the covariance terms. Consequently, the accuracy of $9, in 
(C1) depends on the covariance magnitudes. For example, if 9 
is random, then (C1) would not converge rapidly. Holstrom 
(1963), however, has observed considerable redundancy in the 
atmosphere for geopotential. 

In Obukhov's (1960) method, the eigenfunctions of the 


integral operator, 


p 
{2 B(p.p')o(p')dp' = upa(p), (C4) 
= 


produce the optimum choice for basis functions, o,(P), which 
produce least mean square error for all values of n in (Cl). 
The eigenvalues of (C4), up, measure the variance that their 
associated eigenfunctions represent. Therefore, the order- 
ing of Eneer Uncen TONS 15 Made so that w is in descending 
order. Also, the averaged normalized variance of % that is 


represented by $ is simply 


(C5) 


=) 


y 
~ 


~~ ~ 
um sgpums 
mad 


— 


To extend Obukhov's (1960) procedure to finite differ- 
ences, the independent function is represented by 


bi og.k = PC 14x, jay, kap), (C6) 





where Ax, Ay, and Ap are the grid spacings and i,j,k the 
grid location. Then the autocovariance becomes 
M N 


B =e) » 


ap/MN (C7) 
oe Ey iS] 


nee i adiGm 


The integral operator (C4) now has the form 


m - UL®> 3 


where 2 and m are vertical indexes. The optimum basis 
functions are the eigenvectors of Bom arranged so that 
corresponding eigenvalues are in descending order. There 
are K modes or eigenvectors in this system, which correspond 
to the size of the square array Be,. However, the 
motivaticn of this approach is to allow a partial sum to be 
used that maintains most of the accuracy of the original 


mere tion %1,5,p° 
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